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The calculation of imaginary time displaced correlation 
functions with the auxiliary field projector quantum Monte- 
Carlo algorithm provides valuable insight (such as spin and 
charge gaps) in the model under consideration. One of the 
authors and M. Imada p| have proposed a numerically stable 
method to compute those quantities. Although precise this 
method is expensive in CPU time. Here, we present an al- 
ternative approach which is an order of magnitude quicker, 
just as precise, and very simple to implement. The method is 
based on the observation that for a given auxiliary field the 
equal time Green function matrix, G, is a projector: G 2 — G. 
PACS numbers: 71.27.+a, 71.10.-w, 71.10.Fd 
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where s denotes a vector of HS fields. For a Hubbard 
interaction, one can for example use various forms of 
Hirsch's discrete HS decomposition [j]||. For interac- 
tions taking the form of a perfect square, decompositions 
presented in || are useful. 

The imaginary time propagation may now be written 
as: 



e -2SH = J2 f/-(29, 0) + 0((At) 2 ) where 
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For a given Hamiltonian H = J2 x , y c t T x : yC y + Hi and Ug(2Q, 0) = TT 



,-ATfl,/2 E« „ c t D *,v(.*»)cv -ArHt/2 



its ground state |^o)i our aim is to calculate 



Gx,y( T ) 



{^ \cI(t)c x \^o) 



r >0. 
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Here, c\ creates an electron with quantum numbers x, 



c x (t) = e-( H -^c x e 
tial n = E*+<- 



E Q . 



t(h fiN) ^ anc j ^e chemical poten- 
Hi corresponds to the interaction. 
Within the projector Quantum Monte Carlo (PQMC) al- 
gorithm, this quantity is obtained by propagating a trial 
wave \^t) function along the imaginary time axis [0-0: 



(*o\4(t)c x \* ) 
<*o|*o> 
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(y T \e- 2eH \y T ) 
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lim G<(e,e + r). 
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The above is valid provided that: (^oI^t) / 0. 

To fix the notation, we will briefly summarize the es- 
sential steps required for the calculation of the right hand 
side (rhs) of the above equation at fixed values of the pro- 
jection parameter O. A detailed review may be found in 
p] . The formalism - without numerical stabilization - to 
compute time displaced correlation functions follows Ref. 
P . The first step is to carry out a Trotter decomposition 
of the imaginary time propagation: 



-20H 



e -ATfl 1 /2 e -ArH, e -ATH,/2 



0((Ar) 2 ). 



Here, H t (Hi) denotes the kinetic (interaction) term 
of the model and to At = 20. Having isolated the 
interaction term, Hi, one may carry out a Hubbard 
Stratonovitch (HS) transformation to obtain: 



n=l 



The HS field has acquired an additional imaginary time 
index since we need independent fields for each time in- 
crement. 

The trial wave function is required to be a Slater de- 
terminant: 



N„ 



W = I1 E C ^,n io>. 
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Here N p denotes the number of particles and P is an N s x 
N p rectangular matrix where N s is the number of single 
particle states. Since Ug(2Q, 0) describes the propagation 
of non-interacting electrons in an external HS field, one 
may integrate out the fermionic degrees of freedom to 
obtain: 

g< (e, e + T ) = J2 w sGf (e, e)B ? (e, e + T ) (i) 



where we have omitted the (At) 2 systematic error pro- 
duced by the Trotter decomposition. In the above equa- 
tion, 

S s -(9 2 ,0i) = 

n7=n 1+ i^ ATT/ ^ D{s " n) e- ATT/2 if e 2 >ex 
Bz 1 (& 1 ,e 2 ) if e!>e 2 



(3) where tt-iAt = 9i and u 2 At = 2 



M s =P T B s (2Q,ti)P, W s -. 



det(M s -) 
£ a -det(M ? ) 



and 



G<(e,e) = i? ? (e)[£ ? (0)^(e)]- 1 i ? (0), 

Rg{&) = Bg(&,0)P, L S {Q) = P T B g {2Q 1 Q) 

Restricting ourselves to models where W s is positive 
definite (such as the half-filled Hubbard, half-filled Kondo 
lattice or attractive Hubbard models) we can sample the 
probability distribution with Monte Carlo methods. For 
each auxiliary field configuration we then have to evalu- 
ate the quantity G$(0, Q)Bg(Q, + r) in a numerically 
stable and efficient way. This corresponds to the subject 
of the paper. 

At first glance it is clear that the evaluation of 
G5(0, 0)5^(0,0 + r) is a numerically ill posed prob- 
lem. We illustrate this by considering free electrons on a 
two-dimensional square lattice. 



H=-t V clc T . 



(8) 



<i,j> 



Here, the sum runs over nearest-neighbors. For this 
Hamiltonian one has: 

(*ol4Wcjfl*o> = (*ol4 c fel VI/ o}exp (r(c £ - //)) , (9) 

where e^ = —2t(cos(ka x ) + cos(ka y )), a x , a y being 
the lattice constants. We will assume |\l/o) to be non- 
degenerate. In a numerical calculation the eigenvalues 
and eigenvectors of the above Hamiltonian will be known 
up to machine precision, e. In the case eg — /i > 0, 

(^olcjCgl^/o) = 0. However, on a finite precision ma- 
chine the later quantity will take a value of the order of 
e. When calculating (^o|cjj(T)cd\I) r o) this roundoff error 
will be blown up exponentially and the result for large 
values of r will be unreliable. 

In the PQMC approach and since for a given HS con- 
figuration, we have independent electrons in an exter- 
nal field a similar form is obtained for the time dis- 
placed Green function. The Bg matrix plays the role 
of the exponential factors, and contains exponentially 
large and small scales whereas Gf (0, 0) contains scales 
bounded by order unity. Since we equally expect the re- 
sult G$ (0, + t) to be bounded by order unity, we will 
eventually run into numerical problems when r becomes 
large. 

In order to circumvent this problem, Assaad and Imada 
PI proposed to do the calculation at finite temperatures 
and then take the limit to vanishingly small tempera- 
tures. For the example of free electrons this amounts in 
doing calculation via: 



(*o|4( T ) c fcl*o) = Jim 



exp (r( £ g - jj)) 
fiZZo i _|_ exp (/3(eg - (j,)) 



(10) 



Even if the eigenvalues are known only up to machine 
precision, the right hand side of the above equation for 
large but finite values of (3 is a numerically stable oper- 
ation. To implement this idea in the QMC method, As- 
saad and Imada considered a single particle Hamilton Hq 



which has the trial wave function, \^t) as non-degenerate 
ground state and then compute: 



G5(0,0 + r) 



(11) 
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Tr( e -^^(20,0)ct(r) C:c ^(0,O)) 
T)r(e-/ M *> 17,(20,0)) 



Although the rhs of the above equation may be com- 
puted in a numerically stable way, the approach is cum- 
bersome and numerically expensive. In particular, for 
each measurement, all quantities have to be computed 
from scratch since the ad-hoc inverse temperature j3 has 
to be taken into account. 

Here, we propose an alternative method. We will 
again start with the example of free electrons. Since, 
{^ \cUt)c^\^!o) = 1, 0, we can rewrite Eq. || as: 

(*ol4W c ifl*o) = (<*ol4csl*o) exp ((eg -/*)))' (12) 

which involves only well defined numerical manipulations 
even in the large r limit. 

The implementation of this idea in the QMC algorithm 
is as follows. First, one has to notice that the Green 
function G$(0,0) is a projector: 



G5(0,0) 2 = G5(0,0). 



(13) 



This is simply shown by carrying out a singular value 
decomposition of the N s x N p Rg- (0) and Lg (0) matrices 

Rg(Q) = U rv gD ry gV ry g 

Lg{Q) =VifD lt3 Ui,r. 

Here U r (Ui) is a N s x N p (N p x N s ) and column (row) 
orthogonal matrix, Di r are diagonal N p x N p matrices 
and Vi_ r unit upper triangular N p x A^ matrices. For the 
equal time Green function only U r g, Ui g are important 

G<(0,0) = U^gp^gUrs}- 1 17,,* 
Eq. O then follows from: 

(G5(0, 0)) 2 = UrJp^gUrA' 1 Pl,sU r s] [Ul^r,^ Ul,g 

= G5(0,0). 

This in turn implies that Gf (0i , 03) obeys a simple com- 
position identity 

G<(0 1 ,0 2 )G<(0 2! 3 ) = G<(0i,03) (14) 

since 

G5(0 1 ,0 3 ) = G5(0 1 ,0 1 ) Bg(Q 1 ,0 3 ) = 

(G5(0 1 ,0 1 )) 2 B s -(0i,0 3 ) =G s <(0i,0i)G s <(0i,0 3 ) 
= G5(0 1 ,0 2 )G<(02,0 3 ). 



Using this composition property ( |14| ) we can break up 
a large r interval into a set of smaller intervals of length 
t = Nt\ so that 



JV-1 

G s < (9, 6 + r) = Y[ G< (9 + [n + 1] r 1; 9 + nr x ) (15) 

n=0 

The above equation is the generalization of Eq. 03. If n 
is small enough each Green function in the above product 
is accurate and has matrix elements bounded by order 
unity. The matrix multiplication is then numerically well 
defined. 

We illustrate the efficiency of the method for the 
Kondo lattice model: 



Hklm = — t 



E 

{i,l),o 



C-. C~r 
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Here i runs over the L 2 -sites of a square lattice, (i,j) 

corresponds to nearest neighbors, cl creates a Con- 
ner 

duction electron with z-component of spin a on site i 
and periodic boundary conditions are imposed. Si = 
(1/2) JT , /! fi a , a ih , with a the Pauli matrices. An 
equivalent form holds for the conduction electrons. A 
constraint of one fermion per /-site is enforced. As shown 
in Ref. |U| at half-filling, the PQMC method may be used 
to carry out sign- free simulations of the model. 
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FIG. 1. Imaginary time displaced on-site spin-spin (a) and 
Green function (b) correlation function. We consider a 6 x 6 
lattice at half-filling and J/t = 1.2. In both (a) and (b) results 
obtained form Eq. (|]J) (A) and ( |il| ) (v) are plotted. 

Fig. ffl plots the on-site time displaced spin-spin corre- 
lation functions as well as the on-site Green function for 
a 6 x 6 lattice at J/t =1.2 and half-band filling. Here, 
we consider the total spin: S^ = Si + SS . Both methods 
based on Eq. ( |l5| ) and Eq. (^lj) produce identical results 



within the error-bars. (Had we used the same series of 
random numbers, we would have obtained exactly the 
same values up to roundoff errors which are of the order 

io- 8 ) 

The important point however, is that the method 
based on Eq. (|15|), for this special case, more than an 
order of magnitude quicker in CPU time than the cal- 
culation based on Eq. (11). A calculation following Eq. 
( |l5| ) involves matrix inversions (multiplications) of size 
N p x N p ( (N p x N) ■ (N x TV)). Here, N denotes the 
number of sites. To this we have add that many quan- 
tities required for the calculation are at hand during the 
simulation and do not have to be recalculated. Gn the 
other hand, the method based on Eq. (yd]) involves ma- 
trix inversions and multiplications of size up to 2jV x 2N 
M. In this approach and for given set of HS fields all 
quantities have to be computed from scratch. 

In summary, we have described an efficient method 
for the calculation of imaginary time displaced correla- 
tion functions in the framework of the PQMC algorithm. 
The method is elegant and easy to implement in a stan- 
dard PQMC code and is an order of magnitude quicker 
than previously used methods. We have demonstrated 
the efficiency of the method in the special case of the 
two dimensional Kondo lattice model. Given the ability 
of calculating efficiently time displaced correlation func- 
tions at arbitrarily large imaginary times enables us to 
pin down charge and spin gaps |10|] as well as quasipar- 
ticle weights pd| . Dynamical propertied may equally be 
obtained after continuation to real time via the Maxi- 
mum Entropy method |12| . 
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